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Abstract. This paper is a review of the dynamics of a system of planets. It includes 
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and Poincare's reduction of them to 3N degrees of freedom with a detailed discussion of 
the non-osculating "canonical" heliocentric Keplerian elements that should be used with 
Poincare relative canonical variables. It also includes Beauge's approximation to expand 
the disturbing function in the exoplanetary case where masses and eccentricities are large. 
The paper is concluded with a discussion of systems captured into resonance and their 
evolution to symmetric and asymmetric stationary solutions with apsidal corotation. 
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1 Introduction 



The discovery of the first extra-solar planet in orbit around a main- sequence star was announced 
in 1995. Since then, the number of known extra-solar planets did not cease to grow. As the ob- 
servations are accumulating, planetary systems with 2 and 3 planets are being discovered. More 
than ten are, presently, known^. As the discoveries are recent and many of the discovered planets 
are at the edge of observational capabilities, the uncertainties on their orbital elements and masses 
are large. It is worth recalling that one of the 2-planet systems previously announced, HD 83443, 
vanished from the lists after new observations failed to show the radial velocity variations pre- 
viously identified with a second planet. Another important example of the current uncertainties 
is the "jump" suffered by the determined eccentricity of HD 82943b. During long time, it was 
listed as ~ 0.4, while a new determination using observations over a long span of time gives only 
0.18. By the same occasion, the mass of HD 82943c became twice bigger than believed before 
(Mayor et al., 2004). These discrepancies should be enough to show us how hazardous is the 
task of getting conclusions from the present data and that we should avoid conclusions critically 
depending on the available data. 

In the current state of art, we are just capable of discovering big planets with not too large 
periods. Therefore, the planets so far discovered are big and most of them have orbits close to 
the central stars. Another characteristic is the large eccentricities of many of them. Even if large 
eccentricities favors discovery, this characteristic is not only due to observational bias and needs 
an explanation. (See Ferryman 2000 for a review of the existing hypotheses.) Large eccentricities 
are considered as the result of early migration processes. It is generally believed that the planets 
did not form at their present observed locations, but were driven by a migration process due to 
tidal interaction of the planets with the discs where they were formed (see Papaloizou, 2003). 
Whether this orbital drift is still at work or not is a matter of debate, although it is more plausible 
to assume that it stopped after the end of the planetary formation stage. These early processes 
were also responsible for having driven the (surviving) systems to very stable conditions in which 
orbit periapses appear close to alignment or anti-alignment. This condition is observed in several 
systems. 

Periapses alignment (or anti-alignment) may occur in resonant and non-resonant systems alike. 

In resonant systems, they are the natural states after the system is trapped into a mean-motion 
resonance (see section 6). At variance, in non-resonant system, they are a consequence of the 
angular momentum variations during resonance crossings without capture (Ferraz-Mello et al. 
in preparation). However, and independently on how they reached this condition, an important 
consequence of this type of configuration is that they constitute a stabilizing mechanism for 
planetary orbits, especially if they have large eccentricities. 

Four extra-solar systems seem to satisfy the resonance condition: Gliese 876, HD 82943, 55 
Cnc and 47 UMa. The first two have planets with periods in a 2/1 commensurability, the third 
in a 3/1, and the later close to a 7/3. With regard to Gliese 876, numerical simulations (Laughlin 
and Chambers 2001, Lee and Peale 2002) seem to indicate that these bodies are actually deeply 
trapped in an apsidal corotation (see section 6.1): They exhibit a libration of both resonant angles 
ai = 2X2 — Xi—vui, and also an alignment of their major axes. Apsidal corotation seems to be the 
natural issue of a capture in resonance in the case of two planets with initially low eccentricities 
(Ferraz-Mello et al. 2003; c/. this paper, section 6). The alignment (or ant i- alignment) of 
periapses has not yet been confirmed in the case of the other planetary systems above mentioned. 

The most conspicuous non-resonant system showing nearly aligned periapses is v Andromedae. 
This system has been the object of many numerical and analytical studies (for references, see 
Michtchenko and Malhotra, 2004). The orbit of the planets c and d in this system are such that 
the distance between their periapses oscillates about zero with half-amplitude ^ 60 degrees and 
period ~ 7260 years. 

^For an up-to-date list see the web page "Extra-Solar Planets Encyclopaedia", by J.Schneider at 
www.obspm.fr/planets and links therein. 
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Figure 1: Main systems of coordinates: heliocentric (left), barycentric (center) and Jacobi's 
(right). 

2 Hamiltonian Equations of the N-Planet Problem 

This section considers the Hamiltonian formulation of the problem of N planets orbiting a star in 
an arbitrary configuration. This is a well-known problem in Celestial Mechanics. However, the 
vast majority of papers in Celestial Mechanics deal with the so-called restricted 3-body problems 
in which only 2 bodies have finite masses. Therefore, some basic topics of the general problem 
need to be remembered. 



Barycentric Hamiltonian Equations 

The barycentric Hamiltonian equations of the N+1 body problem are obtained using the basic 
principles of Mechanics. Let mi {i — 0,1, •••,7V) be their masses. If we denote as Xi the 
position vectors of the N-l-1 bodies with respect to an inertial system, and 11^ = milLi their 
linear momenta, these variables are canonical and the Hamiltonian of the system is nothing but 
the sum of their kinetic and potential energies: 

TV 2 N N 

i/ = r + t/ = ivni_GV V ^ (1) 

i=0 1=0 j=i+l -' 

where G is the constant of gravitation and Ay — |Xi — Xj|. This system has, however, 3{N + 1) 
degrees of freedom, that is, 6 equations more than the usual Laplace-Lagrange formulation of the 
heliocentric equations of motion. The system can be reduced to 3N degrees of freedom through the 
convenient use of the trivial conservation laws concerning the inertial motion of the barycenter. 
There are two sets of variables used to reduce to 3A^ the number of degrees of freedom of the 
above system. The most popular reduction, due to Jacobi, is widely used in the study of the 
general three-body problem and of planetary and stellar systems. A less popular reduction is 
due to Poincare; it was first published in 1897, but Poincare himself did not use it because of 
difhculties related with the definition of the associated Keplerian elements (see Poincare, 1905; see 
next section). It appeared in the literature from times to times and started being more frequently 
used around the eighties (Yuasa and Hori, 1979; Hori, 1985; Laskar, 1990). Hagihara (1970) says 
that it was discovered by Cauchy. 



2.1 Poincare's Reduction to 3A^ degrees of freedom 

In Poincare's reduction, the variables are the components of the heliocentric position vectors 
Xi — Xo and the momenta are the same linear momenta 11^ of the barycentric formulation. 
Hence, 

r, =X, -Xo, P^^^,, (2) 
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(i = 1, 2, • • • , N). The given system has N+1 bodies and we thus need to introduce one more pair 
of (vector) variables. Let them be 



N 



ro = Xo, po = ^n,. (3) 



i=0 

A trivial calculation shows that the variables ri,pj (i = 0, 1, • • • , iV) are canonical. Let us, now, 
write the Hamiltonian in terms of the new variables. The transformations of T and U give, 
respectively, 

2 ^ mi 2 ^ mo 2 mo ^ mo ^ mo 



and 



JV JV JV 



i=l i=l i=j+l ^-^ 

where pi = |pi| and n = |ri| = |Aoi|. 

The reduction of the system is immediate. We note, beforehand, that the variable ro is 
ignorablc. Consequently, po is a constant that, by construction, we set equal to zero. The 
resulting equations may be separated into two parts: 

A. The first pair of equations, corresponding to the subscript 0, is: 

Po = ro = gradpg.H'. (6) 

We note that the second of eqns. (6) gives 



P 



JV 

E— • (7) 



mo mo 



B. The canonical equations in the variables r^, p^, (i 0) are given by the reduced Hamiltonian 

H = H-^-^+y^^. (8) 

9 mn ^-^ m.n ^ ' 



This subsystem has 3N degrees of freedom and is separated from the previous one, since po 
is constant. (We did assume po =0.) 

The Hamiltonian of the reduced system is H = Ho + Hi where 

2 A 



N /-, „2 



i=l 
i=JV j=N 

(-- 

,=1 ,.=,+1 X A,, mo 

and ^ ^ 

Hi = G{mo + mi) A = —. (11) 

mo + rrii 

We note that -ffo is of the order of the planetary masses rrii while Hi is of order two with respect 
to these masses. Then Hq may be seen as the new expression for the undisturbed energy while 
Hi is the potential energy of the interaction between the planets. 



4 



It is worth noting that each term 

is the Hamiltonian of a two-body problem in which the mass-point rrii is moving around the mass 
point mo- In fact, from the Hamiltonian given by eqn. (12), it is easy to obtain the second-order 
differential equation 

Tj = -/ij-p- = -G(mo + mi)-p-. (13) 
One of the canonical equations spanned by Fi is 

'-1^ (14) 

This equation apparently contradicts the statements done after which is the heliocentric radius 
vector and Pi is the barycentric linear momentum. However, it only means that the variation of 
Vi in the reference Keplerian motion is not the actual relative velocity of the t"^ body but 'Pi/Pi- 
This means that, at variance with other formulations, the Keplerian motions defined by cqns. (12) 
are not tangent to the actual motions. To distinguish them from "heliocentric osculating" , when 
necessary, we will use the word "heliocentric canonical" . 

2.2 Action-angle variables. Delaunay elements 

The solution of Hq is a set of N Keplerian motions whose generic Hamiltonian is Fj. The 
purpose of this and the forthcoming section is to obtain the Keplerian elements and the Delaunay 
variables corresponding to the relative coordinates introduced before, which must be used when 
a canonical perturbations theory is constructed using Hq as "unperturbed" approximation. For 
that sake, we have to solve the corresponding Hamilton- Jacobi equation and construct the action- 
angle variables of the given problem. We only give here the more important steps characterizing 
the variables appearing in the definitions of their action-angle variables and in the associated 
Delaunay elements. To do it, the study of the planar case is enough and preferable since the 
rotations necessary when the spatial case is considered, although trivial, introduce many new 
equations. All conceptual questions appear in the planar case and have the advantage of making 
the calculations much easier and thus allow the crucial points to be clearly identified. Once 
the conceptual problems are solved in the planar case, the usual three-dimensional equations 
can be easily adapted to give the elements we are looking for. In the plane, the Hamiltonian 
is separable in polar coordinates. To introduce these variables, let us remember that, in the 
reference Keplerian motion, p = /3r. (For the sake of simplicity, we omit the subscript i in the 
forthcoming equations.) Then 

p = /3 (^fa ri/jbj (15) 

where a, b are the right-handed set of unit vectors at r in the positive directions of the increments 
of r, r, ip are the time derivatives of r, ijj in the reference Keplerian motion. The kinetic energy 
term is, then, 

T=^{r'+r^') (16) 
or, introducing the momenta Pr = ^ and = we obtain 

(17) 

The potential energy term is given by 

Uir) = -f (18) 




5 



and the resulting Hamilton- Jacobi equation is the classical one of the two-body problem with /3 
instead of m and ji instead of G{M + m): 

The solution of this equation is well known and does not need to be reproduced here with all 
details. This equation is separable into: 



Pr = \lmE+f)-^ (20) 

= C. (21) 

C, E are integration constants {E = F is the "energy" and C = r x p is the "angular momentum" ; 
the quotation marks are necessary because of the particular definitions of r and p in the considered 
formulation). 

The actions associated with the given Hamiltonian are 

Jr = ^ ^ Prdr J^ = ^ p^dip (22) 

whose integrations give 



Jr = -C + i^PJ^ J^=C. (23) 



-2E 

The Delaunay actions are: 

L = Jj. + = Py/JIa 



G = = /Sy^/ZavT— 

where a and e are two constants introduced in the integration giving the action J^: 
• The mean distance (or semi-major axis) 

def M/? 



(24) 



2E ^''^ 



• The eccentricity 



def L^'^EC^ 



Since, in general, the planets do not move in the same plane, we have to introduce the 
inclinations / of their planes of motion over a fixed reference plane and add the third Delaunay 
action H = l3^/]Ia\/l — cos /. The Delaunay angles £,uj = ^ — v (and fl) are obtained in the 
usual way. 

2.3 Canonical heliocentric elements 

For each planet, we may transform r, p into the elements a,e,X,zu using the same trans- 
formations used to define the ordinary osculating heliocentric elements aoso Soso Aqso'a'osc of the 
two-body problem as functions of m,G{M + m),r, mr. However, the equations giving the os- 
culating heliocentric elements depend on m only through /i. In order to use always the same 
routines, the above equations may be transformed. We substitute, in eqns. (25) and (26), E and 
C by their definitions E = F and C = r x p. We obtain the well-known equations 

a = . . (27) 
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(28) 

where we have used the velocity in the reference Keplerian motion 

w=|, (29) 

instead of the actual planetary velocity, and w = |w|. 

The Keplerian motion corresponding to Hq in Poincare's relative canonical coordinates may 
be obtained with the ordinary routines substituting the heliocentric velocities by 

TTl 

w = -V (30) 

where V is the absolute (i.e. barycentric) velocity. 

The angles are obtained with usual equations. In the planar problem, the true longitude {(p) 
is given by the angle formed by the radius vector with the first axis of the reference system (to be 
obtained through arctan y/x where x,y are the components of r). In the spatial problem, some 
rotations are necessary beforehand. The anomalies may also be easily obtained, starting with the 
eccentric anomaly (u), which is given by 

/ fa r.w \ , 

u = arctan . / . (31) 

Ha-rJ 

The true (v) and mean (£) anomalies are obtained by means of classical 2-body equations. The 
other angles to determine are the longitude of periapsis {u = <j) — v) and the mean longitude 

{x = e + Lo). 

The elements of the reference Keplerian orbit at the time t are a, e, oj, X. Since the parameter 
A is variable, it is convenient to substitute it by the so-called "mean longitude at the epoch" (Ao), 
which is the value of A at a standard "epoch" to : 

X = Xo + n{t- to) (32) 

where n = \ is the mean-motion in the reference orbit. 

2.4 The Conservation of the Angular Momentum 

If the only forces acting on the N-|-l bodies are their point-mass gravitational attractions, the 
angular momentum is conserved: 

N 

£ = ^ mjXi X Xj (33) 
Since irii'X.i = J2o f^i^i = 0) the above equation gives 

N 

£ = ^riXpi, (34) 

that is 

N 



C = ^l3iJfi,at{l- ef) -ki (35) 



where kj are the unit vectors normal to the Keplerian planes. This is an exact conservation law. 
In this equation Ui and Ci are not the usual heliocentric osculating elements but the canonical he- 
liocentric elements defined by equations (27) - (28) where Wj are the absolute velocities corrected 
by the factors mj//3j. 
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The conservation law given by eqn. (34) is also true if Jacobian coordinates are used. However, 
the expression 

£ = ^TOi^/ijaj(l - e|).kj (36) 

where ai and Cj are the heliocentric osculating elements (Keplerian elements defined by eqns. (27) 
- (28) with the heliocentric velocities instead of Wj) often found in the literature, is not a true 
conservation law. One may easily see that: 

N 

£ = £ - ^ m^Xo X Xo (37) 

i=l 

showing that the quantity C has in fact a variation of order 0{m^). 



2.5 Two-body Expansions 

For the sake of future calculations, let us recall some series expansions of the two-body problem. 
These expansions are helpful in the task of writing computer codes for automatic expansion of 
Hi and hold in all systems of elements founded on unperturbed Keplerian motions. 

The first result to be recalled concerns the Fourier expansion of some functions of the radius 
vector and true anomaly. They are the convergent series 

\ n oo 

- cos(fc/) = ^(x;''=+X!!''=)cos(i£) (38) 



j=0 



sin(fc/) = J2i^^''-X^f)sm{j£) 

3=0 



where the superscript n may be either positive or negative. The coefRcients X^'^ are the Hansen 
coefficients (see Tisserand, 1960; Kaula, 1962). Hansen coefficients are functions of the eccentric- 
ity. They may be developed into power series of the eccentricities: 

oo 

^;'' = e"=-^'E^X>«+«.«'^ (39) 

=0 



(wi = max(0, j — k) and U2 = max(0. A: — j)) where the numbers i^J+Jj^ are the Newcomb 
operators. Newcomb operators obey to some simple recurrence relations, which allow them to be 
easily calculated for all values of the indices (sec Brouwcr & Clemence, 1961). 
Introducing eqn. (39) into eqn. (38), we obtain, after some algebra. 



)n oo oo 

C0s(fc/) = X! X! ^n,k,i,me' COS (mi) (40) 
i=0 m=— oo 



oo oo 

r 



sin(fc/) = X] X] Cn,k,i,me' sin (ml) 



i=0 m=—oo 



,k,i.ra are constant coefficients expressed as functions of Newcomb operators. 
These coefficients, first calculated by Leverrier, do not depend on the orbital parameters and may 
be calculated once for all. They have some interesting properties. The most important of them 
is the d'Alembert property: Bn,k,i,m = Cn,k,i,m = when |m| < z or when |m| — i is odd. 

The latest expansions are power series in e and their convergence depend on the singularities 
of the analytic function u = u{e,£), which are at |e| = 0.6627434- • •. This is the convergence 
radius of the given series (see Wintner, 1941). 
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3 Expansion of the Disturbing Function 



The Hamiltonian equations in relative coordinates may be used to study the planetary motions. In 

analytical studies, once introduced the new variables, the next step is to write Hi in terms of the 
Keplerian elements. A well-known approach to this problem is the classical Laplacian expansion of 
Hi into a Fourier series in the angles and a power series in the eccentricities, which introduces the 
fimctions of the semi- major axes known as "Laplace coefficients" . Another expansion sometimes 
found in the literature uses the expansion of in Legendre polynomials of the ratio of the 
distances of the two planets to the central star. These expansions work well in their domains of 
validity. The Laplacian expansion is a good approximation if the orbital eccentricities arc small. 
However, the radius of convergence of the expansion decreases (see Ferraz-Mello, 1994) with the 
increase of the ratio a of the two semi-major axes. For a ~ 0.6 the series is no longer convergent 
for eccentricities as small as ^ 0.2. The expansions with Legendre polynomials arc more stringent: 
they may only be used in the study of well hierarchized systems where the ratio of the distances 
of the perturbed and perturbing bodies to the central body remain small forever. This is the case 
of the lunar theory, in which the motion of the Moon around the Earth is disturbed by the Sun. 
Otherwise, the convergence of the expansion in Legendre polynomials is very slow and its use in 
planetary problems accounts for many wrong results. We present, in this lecture, an improvement 
of the technique first developed by Beaugc (1996). This "expansion" is valid in large domains of 
the phase space excluding a domain around the singularities associated to collisions between the 
two bodies. In Beauge's approximation, the number of terms necessary to represent Hi depends 
on the magnitude of Hi in the domain to be studied: Near the minimum of \Hi\, a few terms are 
enough to have a good representation. This number increases quickly as we approach orbits that 
may come close to a collision. At variance with Beauge's early expansion, the present one (Beauge 
& Michtchenko, 2003) has no explicit restrictions with regard to eccentricities and inclinations. 

3.1 Beauge's approximation. The parameter S 

The big problem in the expansion of Hi comes from the terms having A in denominator. In 
heliocentric coordinates, we can write: 

1 = (r2 + ri-2rir2COs5)-'/' (41) 

where S is the angle between both bodies as seen from the central mass. Introducing the ratio 
p = ri/r2, eqn. (41) becomes 

^ = (l + p2_2pcos5)"^^^ (42) 

Instead of expanding this function in Fourier series of S (Laplace approach) or power series of p 
(Legendre polynomials), a best-fit approach is used. We write 

where 

x^p^-2pcosS (44) 
and represent the function (1 + a;)~^/^ by a polynomial of order in x: 

N 

(l + a;)-i/2^^&„x" (45) 

n=0 

whose coefficients fo„ are determined numerically through a linear regression. 

The variable a; is a measure of the proximity of the initial condition to the singularity in . 
It is equal to —1 at the singularity, and takes values larger than this for every point (p, cosS') 
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Figure 2; Relative error of the approximation of l/vT+x given by eqn. (45) in function of x, for 
two values of 5. The continuous line shows the case 5 = 0.1 and the broken line shows 6 = 0.01. 
In both examples, A'' = 30. 

away from the collision curve (see fig. 3). Wc note that the values of p and S are not separately 
significant; only the distance from the singularity is important. 

The numerical fit is performed using values of .t > ~l + S, where 5 is a positive parameter 
close to zero. The smaller its value, the better the approximation to the real function near 
the singularity. However, when 6 is small, the number N of terms to be considered in the 
representation of (1 + .t)^^/^ to guarantee an adequate precision for all values of the independent 
variable is necessarily large. 

Figure 2 shows the relative error of (45) for N = 30 and two values of 5. We can see that for 
most of the interval of .x, the fit with (5 = 0.1 yields a miich higher precision than the fit with 
S = 0.01. In the fit with 5 = 0.1, the errors are of the order of 10~^, that is, about 3 orders of 
magnitude lower than in the other case. Conversely, as a; ^ —1, the fit with 5 = 0.01 is more 
precise. Larger values of N will diminish the error in both cases, but at the cost of increasing the 
number of terms enormously. 

In the general case, the motion of the two bodies is unconstrained and the distance between the 
two planets is minimum in a symmetric conjunction with the outer planet at the periapsis and the 
inner planet at apoapsis. In this case, we have to choose S < (1 — ^)^ where ^ = a(l+ei)(l — 62)"^. 
Beauge's technique no longer requires that the ratio of the distances of the two planets is small, 
but it requires ^ < 1. However, when the planets are in resonant motion, the method is valid even 
for crossing orbits because the resonance does not allow the planets to come close one to another. 
The limits of x when the motion of the two planets is constrained by a 2:1 commensurability 
(a = 0.63) are shown in figure 3 in the particular case where 62 = 0. 

The geometry of the curves in fig. 3 follows very closely (but not identically) the topology 
of the phase portrait of the 2:1 resonant restricted three-body problem averaged over short- 
period terms. The maximum value of x^i^ lies at ei = 0.8, (Ti = (on the horizontal axis) and 
corresponds to the minimum of |. This point is very close to the corotation stationary solution 
of the 2:1 asteroidal resonance (ei = 0.73 when 62 — > 0; see Ferraz-Mello et al., 1993). Similarly, 
the minimum value of x^nin (equal to —1) corresponds to the singularities of Hi. There is no 
direct relationship between the eccentricity and Xmin ■ An orbit with a large eccentricity near the 
corotation center may have a larger value of Xmin, while an almost circular orbit with a lower 
eccentricity may reach values very close to the limit x = —1. It is worth recalling that several 
extra-solar planet pairs observed in resonant configuration lie near corotation centers where Xmin 
is large and good Beauge's approximations may be obtained with small N. 
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Figure 3: Limits of validity of Beauge's approximations for planets in the 2:1 resonance with 
62 = for different values of 6. The thick black line on the left-hand side is the locus of the 
points where .Tmin = — 1 (collision curve). The non- labeled curves adjacent to it correspond to 
S = 0.001. Horizontal axis: ei coscti; Vertical axis: ei sinai. (cti = 2A2 — Ai — wi). 



3.2 The Direct Part 

To transform the above approximation into a function having the form needed in a theory, many 
transformations have to be done. Introducing the explicit expression for x into eqn. (45), it 
becomes 

g:.^^c.(-2)^- J p^'^-^cos^^ (46) 

fe=0 j=0 \ / 

where the Ck are constant coefficients, easily obtainable in terms of the bk- 

From now on, wc will restrict ourselves to coplanar orbits. Changing from powers of the cosines 
to multiples of the argument, and introducing the planar approximation S = fi — f2 + Azu, we 
can rewrite it as: 

N N — k ^ \m/ \ —m—1 

^:..^5]2^fc,,a'"P p cosk{f,-f2 + Aw) (47) 

where m = 2i + k. 

At last, introducing eqn. (40) into the expression of the direct part of the disturbing function, 
and reordering the terms, we get: 

00 00 Af N — l 

J2 X^X!^'>*-^2i+i,j,fc,m,na2*"^'ele^cos(m£i-n£2 + ^AK7) (48) 

j,k=0 m,n= — 00 1=0 i=0 

where the coefficients D2i+i,j,k,m,n are given by: 

D2i+l,j,k,m,n = {B2i+l,l,j,\m\ + sign{m)C2i+Llj,\m\) X (49) 

{B-2i-l-l,l,k,\n\ + sign(n)C_2i_;_i,;,fe,|n|) 

and 7m is a simple bi-valuated function defined as: 

r 1/2 if m = 

1 if m>0. (^0) 
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Eqn. (48) multiplied by the factor gives the term of the direct part corresponding to 

the given pair of planets. 



3.3 The Indirect Part 

In Poincare heliocentric relative coordinates, the indirect part of Hi is (see eqn. (10)): 

N N 



i=l j=i+l 



Too 



(51) 



The linear momenta pi may be obtained from the derivatives of the vector radii ri{t),Vj{t) in the 
Keplerian reference orbit (see eqn. (14)). Then, 



N N 
1=1 j=i+l 



Vi{t)Vj{t) 



Too 



or 



rr. /3l/32 

Ti = 711712 

TOO 



dx^dx2 dyi dy2 
d£i d£2 d£i d£2 



(52) 



(53) 



where £i are the mean anomalies and Ui the mean motions. Xi and are the components of 
and are given by Xi = Vi cos(/j + vji) and yi = ri sin(/i + Wi). 

In the sequence, we substitute the mean motions by the values issued from Kepler's third law 
and put into evidence the same factor used at the end of the previous section. Hence 



02 



d { xi\ d 



X2 



+ 



d f yA d / 2/2 



d£i \ai J d£2 \a,2 J d£i \aij d£2 \ 02 



(54) 



A = 



j3i/32_ « 1 — jg hereafter equal to 1, introducing an error of third order in 

mim2 2mo ^ ' ° 

the planetary masses. Using the expansions given in section 2.5, there follows 



Xi 
Oi 

yi 

ai 



00 00 



= E E It,je'cos{j£i+wi) 

i—0 j— — oo 



00 00 



E E -fi,je*sin(j£i + OTi) 

i=0 j=—oo 



where 



1 
2^ 



Bi,i.t^\j\ +sign(j)Ci_i,,jj| 



(55) 



(56) 



After the differentiation of these equations with respect to the mean anomalies, and substitution 
in Ti, we obtain 



Ti 



Gmim2 
172" 



02 a 



,7ie|e2 cos {m£\ — n£2 + Akj). 



(57) 



j,k—0 m,n= — 00 



Notice that, except for the dependence on a, this series is formally similar to that giving the 
direct part of F\. To complete the similarity, we can substitute the factor a"^/^ by a power series 
expansion in the neighborhood of the exact resonant value and write it as 



2N 



-1/2 



(58) 



i=0 
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where Ai are constant coefficients. With this change, Ti now reads: 

^ 27V oo oo 

Ti = — - 'AimnIj,mIk,nCt''e{e2 cos {mii - ni2 + Aro), (59) 

(l2 ^ ^ ^ 

i=0 j,fc=0 m,n= — oo 

which is the final expression for the indirect part of the disturbing potential. 
3.4 The Disturbing Function 

Since both parts arc formally similar, wc can unify both expressions and obtain a single series for 
the disturbing function of the planetary three-body problem in heliocentric relative coordinates: 

^ oc ca N 2N 

H, = y y yy R^,j,k,m,n,la'eie^ cos {mh - + 1 Aw) (60) 

do ^ ^ ^ ^ 

j,k=0 m,n=—oo 1=0 i=0 

where a = — and 

Ri,j,k,m,n,l = -^l,{i-l)/2DiJ^k,m,n — Sl,0-AimnIj^rnIk,n', (61) 

Sifl is Kronecker's delta function. Note that these coefficients are constant for all initial conditions, 
and therefore need only be determined once. (For more details, see Beauge & Michtchenko, 2003.) 

It is important to note that each term in Hi depends on the mean anomalies ii and on the 
diflference of the periapses longitudes Aw. This means that if the arguments are written in terms 
of longitudes A,, to, only, they become KiAi + K2X2 + + K4W2 with ^ = 0. That is, Hi 
is invariant to rotations of the reference frame. 



4 Secular Dynamics of 2 Planets 



The study of the secular dynamics is the study of the secular part of the Hamiltonian, obtained 
after an averaging over the mean longitudes. We will restrict ourselves in this text to the case of 
only two planets. To the first-order of the masses, the averaged Hamiltonian is the mean value 
of H: 

1 

<H > = - 

L 

or 



<:H> = ^ I I Hd\id\2 (62) 



2 



<H> = -Y,^- RseciLi, Ii, Aw) 

where we have introduced the Delaunay variables Li,Ii = L^ — Gi defined in section 2.2. Because 
of the invariance of Hi with respect to rotations, once the are averaged o\it. only terms with 
arguments n^wi + K4W2 with ^3 = —K4 can remain in .Rsec- That is, < H > depends on only one 
angle, the difference Aw. This means that the averaged equations have three ignorable angles, 
that is, three first integrals (conservation laws). They are 

Li = const. 

L2 = const. 

K2 — Ii + I2 = const. 

The third of these integrals, 

K2 = Ii+h= Li{l - ^1 - el) + 1.2(1 - ^1 - ei). (63) 

was called Angular Momentum Deficit by Laskar (2000). It is a combination of the conservation 
of the angular momentum (Gi + G2=const., in the planar case) and the secular invariance of 
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Figure 4: Secular variation of the osculating eccentricities of vAnd planets. The eccentricities of 
the two largest planets have variations in anti-phase as necessary to have AMD w const. 

the Lj. The /j are positive quantities increasing from when = to Lj when Cj = 1, that is, 
< 7j < ij. Therefore, in a system formed by only two planets, the /j shall vary in contrary 
directions and so shall vary the eccentricities: When one eccentricity increases, the other decreases 
(see fig. 4). 

This conservation law has some algebraic consequences. Assuming that ai < a2, we have the 
following possibilities: 

• K2< Li < L2 

1. Ii and I2 cannot reach their maximum values Li and L2, resp.; 

2. ei < 1 and 62 < 1 (for all t); 

• Li< L2< K2 

1. The AMD does not bound the eccentricities (both can reach e = 1); 

2. ei > and 62 > (for all t). 

• Li< K2< L2 

1. /i can reach its maximum value Li; 

2. I2 cannot reach its maximum value L2 {I2 < L2 — Li); 

3. The AMD does not bound ei (it can reach ei = 1); 

4. The AMD bounds 62 (62 < 1 for aU t); 

5. I2> K2-Li>0; 

6. 62 > (for all t). 

The conservation law is also found in N-planet systems. In the coplanar case, the angular 
momentum deficit is 



It is worth emphasizing that this conservation law of the averaged system is not a rigorous one 
as the Angular Momentum conservation discussed in section 2.4. It is approximated and valid 



N 
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strictly only as far as the hypotheses done to average the system are satisfied and the semi-major 
axes remain approximately constant. 

The equations of motion derived from < H > are 



h = - 



d<H> 
dAxu 



<i II 

dh 



(64) 



This system has only one degree of freedom and is integrable. 
4.1 The Mode I and Mode II Periodic Motions 

The exact solution of cqns. (64) is not easy to obtain analytically, but some insight can be gained 
by examining their equilibrium points (which correspond to periodic solutions of the two-degrees- 
of- freedom Hamiltonian < H >). They are defined by 



For non-singular Ji (/i ^ and Ji ^ K2), we have two trivial solutions: Aw = and Aw = n. 
These solutions are often referred as Mode I {Aw = 0) and Mode II (Aw = tt). In mode I, the 
lines of apses of the two planets arc aligned having the pcriapscs on the same side; In mode 
II, the situation is similar but the two periapses are in opposite directions (the periapses are 
anti-aligned). Ordinary motions are oscillations around these fixed points. 

Solutions of the above equations corresponding to the masses, semi-major axes and energy 
level of the planets c and d of fAnd are shown in Fig. 5. They are presented in two different 
planes. One in which the coordinates are eicosAsj, eisinAro (ei is the eccentricity of vAnd 
c) and another, not independent, in which the coordinates are e2C0sAro, e2sinAro (62 is the 
eccentricity of vAnd d). 

On each figure, we can see the two fixed points above called Mode I and Mode II. In the 
left-hand phase plane, corresponding to the eccentricity of planet c, motions aroimd the Mode 
I fixed point dominate; the Mode II fixed point lies near the left-hand boundary of the energy 
surface. In the right-hand figure, corresponding to planet d's eccentricity, the situation is reversed 
and the flow is dominated by motions around the Mode II flxed point which lies near the center. 
(For a discussion on the periodic orbits corresponding to the fixed points, see Michtchenko & 
Ferraz-Mello, 2001; Michtchenko and Malhotra, 2004.) 

It is important to note that even though the motion of angle Aw may be cither an oscillation 
(about or 180°) or a circulation, there is no separatrix associated with an unstable infinite- 
period solution separating these motions. To better understand this feature, we plot by dashed 
lines two special solutions on each figure. These solutions arc associated with the singularities 
in eqns. (65), which take place at /i = and Ii = K2 (that is, /2 = 0). One of these solutions, 
presented by the curve S\, was calculated with initial condition 7i ~ and is seen as a smooth 
curve passing through the origin on the left-hand side figure. At variance, ^2, calculated with 
initial condition /2 ~ is seen as the 'false' separatrix between the domains of the motion around 
the two different fixed points. An analogous situation is seen in the right-hand-side figure where, 
now, 5*2 is a smooth curve passing through the origin and Si separates the domains of the motion 
around the two different fixed points. The motion along these separatrix- like curves is such that, 
when the representative point in one plane passes through the origin, in the other plane it is at 
the boundary of the separatrix-like curve and jumps from one boundary to another. However, 
such jump does not mean that the motion is passing through a singularity. It is just the result of 
the topological inadequacy of the plane to represent these solutions; they would be better drawn 
over a sphere (see Pauwels, 1983). 

Figure 5 shows some important features of the secular motion of two planets. In solutions close 
to Mode I (the right-hand side fixed point), the secular angle Aw oscillates about and the planet 
eccentricities undergo small oscillations aboiit the value corresponding to Mode I equilibrium. In 
a similar way, the solutions close to Mode II (the left-hand side fixed point), the secular angle 



/i -0, 



Aw = 0. 



(65) 
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Figure 5: Secular variations of a system of planets with the same masses and semi- major axes 
as planets c and d of vAnd. The outer curves show the boundary of the energy manifold. The 
dashed curves 5*1 and 5*2 represent motions through the singularities of eqns. (65) (see text). The 
actual motion of planets uAnd d and vAnd c (see fig. 4) is an oscillation around Mode I fixed 
point. 



Am oscillates about 180°. At mid-way from Mode I to Mode II, there is a large region of the 
phase space, corresponding to solutions where the motion of the secular angle Anj is a prograde 
circulation. 

The motions around Mode I and Mode II are two opposite stable ways of the planetary system 
to be aligned. In Mode II, the closest approaches between the planets occur when uAnd c is at 
apoapsis and uAnd d at periapsis, simultaneously. This situation can never occur in Mode I. 

Fig. 5 is akin to surfaces of section of the two-degrees-of-freedom system. The curves in each 
plane are defined by initial conditions on the plane plus one condition out of the plane (the other 
eccentricity, or, equivalently, K2), which is adjusted in such a way that all curves correspond to 
the same energy. Therefore, it is not a phase portrait. (Phase portraits of one-degree-of-freedom 
Hamiltonian are sets of trajectories of different energies. See the phase portraits of < > in 
Pauwels, 1983.) This choice makes fig. 5 more suitable for comparison to similar plots obtained 
for 2-planet resonant systems (Michtchenko & Ferraz-Mello, 2001; Callegari Jr. et al., 2004). 



5 Resonant Dynamics 

In the previous section, the Hamiltonian was averaged over the two mean longitudes Ai and A2. 
This procedure is not valid if the two planets have commensurable periods, since, in this case, Ai 
and A2 are no longer independent: 

P + q T2 p + q 

resonance 



P Ti p 

The averaging over the two longitudes will kill all terms depending on the longitudes including 
those depending on the critical combination 

{p + q)X2 -pAi. 

However, these terms play a major role in the dynamics of the two planets and should remain in 
< H >. To preserve them, we define, before the averaging, the following set of planar canonical 
variables: 
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Ai Ji = Xi + s(Ji + h) 

A2 J2=i2-(l + s)(/l+/2) if.f.. 

(1 + .s)A2 - sAi -wi=ai h=Li- Gi 

(1 + s)A2 — sAi - n72 = (72 h ^ L2 — G2 

where s = p/q. The two angular variables cTj are the critical angles. With the angles thus 
introduced, the generic argument mii — n(.2 + l^w of the disturbing function becomes (m — 
l)ai — (n — 1)(J2 + [m(l + s) — n,s](Ai — A2). Note that, because of the invariancc of Hi to 
rotations, the mean longitude only appears through the mean synodic longitude Ai — A2. It is 
easy to see that the "action" conjugate to the missing angle is the total angular momentum 
C = Gi + G2 = {Li — Ii) + {L2 — I2) = Ji + J2- The averaging over the mean longitudes (or over 
the mean synodic longitude) can, now, be done and the critical angles will be preserved inside cti 
and (T2. 

After the averaging, 

< > = - ^ -^JJ- - Rres 
i=l * 

where 

Gmim2 \ ^ , I. r I I, XT 

■^res = C[...]Q; 6^63 COs[m gCTi + n (0-2 — Cri)J 

o-i . . , . . 

The momenta whose conjugate angles no longer appear m. < H > are first integrals (only 2, 
now): 

Ji = const. 

J2 = const. (67) 

Ji + J2 = Gi + G2 is the Angular Momentum, whose conservation in the system before the 
averaging was discussed. It is worth emphasizing the fact that the Lj (i.e. the semi-major axes) 
are no longer invariant. 

The two integrals above may be combined to give 

(1 + s)Li + SL2 = const. (68) 

(Sessin and Ferraz-Mello, 1984). This integral of the resonant dynamics means that ai and 02 

vary in anti-phase. When one of the semi-axis increases, the other necessarily decreases. 
The above variables may also be combined to give: 

^M£) = 7i+/2 = const. + —. (69) 

s 

The AMD also is no longer invariant, but its variation is small and thus limitations of the 
eccentricities similar to that observed in the secular motion (but different) exist. 

5.1 Resonant Stationary Solutions. Apsidal Corotation. 

The averaged system is, now, an irreducible two-degrees-of-frccdom system. An important feature 
of this system is the existence of stationary solutions (Beauge et al., 2003; Ferraz-Mello et al. 
2003; Lee and Peale, 2003). These solutions are defined by the equations 

dli^ ^ d<H> ^ dRres ^ Q da^^ d<H> ^ ^ 

dt dai dcji ' dt dli 

They arc such that It and at are constant (except for the short period terms eliminated by 
the averaging and for contributions of higher orders). Constant /j's mean semi-major axes and 
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eccentricities constant in these solutions; ai and (T2 constant mean that Azu = (J\—<ji is constant, 
that is, the periapses are moving with same velocities so that their mutual separation do not vary. 
This frozen relative state in resonant systems is known as apsidal corotation. 

Equations (70) may be studied separately. The first equation says that the stationary solutions 
lie at the extrema of the function i?res with respect to the variables cTj. These extrema depend 
only on the ratio of the masses of the two planets and on the eccentricities (constants in the 
stationary solution). The factor '^"^^'"^ does not affect the results. 



e1 = 0.02 ; e2 = 0.02 e1 = 0.02 ; e2 = 0.04 




Figure 6: Contour plots of R^cs in the 2/1 resonance for 4 given pairs of eccentricity values. 
Abcissas: ai = 2X2 — Xi — wi; Ordinates: Aw = vd\ — wi. 

Figures 6 arc contour plots of the function i?i.cs for given ei, ei and a (taken at a = 0,1/02 ~ 
0.63, value corresponding to the resonance 2/1). For the sake of an easier interpretation, we used 
the angles a\^Aw = <T2 — <ti, instead of ci, (T2. The extrema seen in these figures may correspond 
to stable stationary solutions or not. < > is a function of 4 variables and only 2 variables 
are considered in these figures. Therefore, what appears as an extremum in this picture is not 
necessarily one extremum in the full phase space. The stable solutions considered in this paper 
are those corresponding to the centers in the white areas. However, one should be aware that they 
are not the only stable stationary solutions in this problem (see Hadjidemetriou and Psychoyos, 
2003). 

The two first plots correspond to low ei (e\ = 0.02). For small 62 {ei = 0.02) the extremum 
corresponding to stable solutions is such that Aw = tt (ui = 0,(72 = tt). In this solution, the 



18 



periapses are anti-aligned. When 62 is larger (e2 = 0.04 in the right-hand plot), the extremum 
seen in the left-hand plot becomes a saddle point and a bifurcation gives rise to two extrenia sym- 
metric with respect to the saddle. These extrema correspond to asymmetric stationary solutions 
where Anj = (72 ~ ci remains constant but with a value not necessarily equal to zero or tt or 
commensurable with tt. The second row of plots correspond to high ei (ei = 0.2). For small 62 
(e2 = 0.01) the extremum corresponding to stable solutions is such that An? = (fxi =172= 0). 
In this solution, the periapses are aligned. When 62 is larger (e2 = 0.05 in the right-hand plot), 
the same phenomenon seen in the first row occurs: the extremum seen for low 62 gives rise to two 
extrema symmetric with respect to the saddle. As in the previous case, these extrema correspond 
to asymmetric stationary solutions. These asymmetric solutions, depending on the eccentricities, 
may be found on a large set of points tii, Accj. 




Figure 7: Asymmetric stationary solutions. Left: |An7| = 84°, ei = 0.286 and 62 = 0.3. Right: 
\Am\ = 104°, ei ^ 0.17 and 62 = 0.38. 

To complete the determination of the stationary solution, we need to solve the remaining 
equation: 

d<H> ORr^s „ 

___ = - (1 + s)n, + ^ = 0. (71) 

At variance with the previous equation, the solutions of this equation depend on the masses. 
However, it is easy to see that it depends almost only on the ratio of the masses of the two planets. 
Indeed, the are constants in the stationary solutions and the commensurability relation at the 

resonance is 

p 

sui — (1 + s)n2 = 0; 

In the remaining part, the derivatives of i?i.cs with respect to change the dependence on the 
masses. We remember that 



(1 + !!^)-^/^yG^,{i /T^). 



mo 

Thus, the coefficient in front of the summation in i?].cs becomes linear in the planet masses after 
the derivative with respect to li. Therefore, eqn. (71) has the form Aiirii + A2m2 = 0, whose 
solutions do not depend on the masses themselves but only on the mass ratio 7712/1111. This is 
not a rigorous statement. In fact, the semi-major axes and eccentricities are functions of li that 
include also the factor mp -I- m^. This means that Ai and A2 are independent on the masses only 
in a first approximation. Even if their variation with the masses is small for the range of masses 
of the considered planets, this variation exists and will affect the solutions in case of large ratios 
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Figure 8: Loci of the stationary corotation solutions of the 2/1 resonance for several mass ratios 
m2/mi. Top figures correspond to the symmetric solutions of the two left-hand side plots in fig. 
6. The points corresponding to two early determinations of the elements of Gliese 876 are shown 
in one of these plots. The bottom figure corresponds to the asymmetric solutions of the two 
right-hand side plots in fig. 6. The line across these curves shows the values of the eccentricities 
for which 0.63(1 -I- ei) = (1 — 62). In all panels, the thick line shows the boundary between the 
domains of symmetric and asymmetric solutions. 



mi/rriQ. Beauge et al. (2003) have shown that the stationary orbits obtained in this section exist 
for planet masses less than ^ 10^^ of the star mass. 

The above equations were used to find apsidal corotation solutions in the case of planets in 2/1 
and 3/1 mean-motion resonances. The relationship between eccentricities and mass ratios in some 
of these solutions are shown in fig. 8. The top panels correspond to symmetric solutions. In the 
left-hand side panel, the periapses are anti-aligned. This is the case of the two innermost Galilean 
satellites of Jupiter: lo and Europa. In the right-hand side panel, the periapses are aligned. This 
is the case of the two planets in orbit around the star Gliese 876. The thick lines in the two top 
panels show the boundary above which symmetric solutions no longer exist. At the thick line, 
the solutions bifurcate into pairs of asymmetric solutions. The relationship between eccentricities 
and mass ratios in the domain of asymmetric solutions is shown in the bottom panel. It is worth 
noting that the mass ratio mi/m2 in the bottom panel is always smaller than a limit close to 
1.0. This situation is often called "exterior case" since it corresponds to have the smaller body 
in an orbit exterior to that of the more massive one. Asymmetric apsidal corotations are known 
in the exterior asteroidal case (Beauge, 1994). Asymmetric periodic solutions in the restricted 
three-body problem were first shown to exist by Message (1958). We may also mention a similar 
behavior, in deep resonance, of the Laplacian critical angle of the Galilean satellites of Jupiter: 
Al - 3A2 + 2A3 (Greenberg, 1987). 
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6 Capture into Resonance 



In this section, we present the results of a series of numerical simulations of the dynamical 

evolution of fictitious pairs of planets under the action of a non-conscrvativc pcrtiubation that 
adds angular momentum and energy to the orbit of the innermost planet. The planets are small 
(some 10~^ of the central body mass) and the mass ratio is m^/mr = 0.538 (i.e., the so-called 
exterior case). The actual calciilations were done with satellites instead of planets, but the 
physical nature of the system does not matter in the following discussion. The physics and used 
methodology are in (Ferraz-Mello et al., 2003). 

The initial distances to the star arc just behind the 2/1 resonance: a = a\ja2 = 0.612. When 
the semi-major axis of mi increases, ai increases and the mean-motion resonance (a = 0.63) 
between mi and m2 is reached. Capture then can take place. The probability of capture depends 
on the rate of variation of ai if the rate is high, the orbit crosses the resonance without 
capture, one phenomenon very well studied in the case of one massless particle. Other factors 
influencing the probability of capture are the orbital eccentricities - capture is more probable 
when orbital eccentricities are small (Dermott et al., 1988; Gomes, 1995). In our calculations, 
initial eccentricities were lower than 0.001 and the physical parameters were adjusted to have 
slow resonance approximation. Figure 9 shows the evolution of the semi-major axes. 



6.1 Capture into Apsidal Corotation 

The system evolves with the innermost orbit receding from the central body (because of the 
non-conservative forces acting on mi) up to the moment where the system is captured into a 
resonance. a2 is almost constant. When the 2/1-resonance is reached, the system is trapped by 
the resonance. As known since Laplace, after the capture, mi continuously transfers one fraction 
of the energy that it is getting from the non-conservative source to m2, so that a2 also increases. 
One may note from fig. 9 that, after the capture into the resonance, Oi increases at a smaller 
pace than before the capture. The increase of the semi-major axes is such that the ratio ai/a2 
remains constant. 






a2 




▲ 



Figure 9: Evolution of the semi-major axes before and after the capture into resonance. Triangles 

mark the moment of the capture. Dashed lines extrapolate the evolution before the capture and 
show the change in slope of the evolution lines, (arbitrary units) 

Figures 10 show the variation of the eccentricities, critical angles (Tj = 2A2 — Ai — tUj and Azu 
in the same time interval as the previous figures. They show that, after capture, the two critical 
angles become trapped in the neighborhood of and tt, respectively and, consequently, the angle 
Atu is trapped in the neighborhood of tt. The capture into a symmetric apsidal corotation with 
anti-aligned periapses is thus simultaneous with the capture into the resonance. 



6.2 Evolution after Capture 

Figure 10 also shows that, after some time, (J2 jumps from tt to and the apsidal corotation 
becomes one with aligned periapses. This change is not the result of a discontinuous process. 
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Figure 10: Variation of the eccentricities, critical angles (ct^) and Aw in the same time interval 
as figure 9. The vertical dotted lines show the moment of the capture. Aw is only shown in the 
final part since it does not differ significantly from a2 in the time interval between the capture 
into resonance and the bifurcation. 



The left-hand side plot shows that the change happens when the eccentricity 62 is zero. Thus, we 
may describe the process by a momentary circularization of the orbit such that, when it becomes 
an ellipse again, the periapses is not at the same side as before. The large transients shown by 
the variation of the angle (T2 are just due to the sensibility of the angle W2 to small changes when 
62 ~ 0. 

The apsidal corotation with aligned periapses does not last long. The figures show that the 
angles depart from zero and the apsidal corotation becomes asymmetric. At this moment, there is 
a discontinuity in the rates of variation of the eccentricities (the elbows seen in the curves ei{t)). 

Figure 11 shows the evolution over a time interval almost 10 times longer. The first point to 
stress here is that such a time span is likely beyond physical signification. Each panel of fig. 11 
combines the variations of ei (resp. 62) and Aw in a same plot in polar coordinates in which the 
radius vector is the eccentricity and the polar angle is Aw. We itemize the important points to 
be noted: 



• The eccentricity ei increases monotonically. When it reaches ^ 0.46, the asymmetric apsidal 
corotation changes back to a symmetric configuration with aligned periapses. 

• After the bifurcation, the followed solution lies in the upper half-plane. This is not the 
only possibility. The motion could have followed a mirror path in the lower half-plane. The 



22 



Figure 11: Polar plots in the planes e, expiAcc. 

probability of following one or another branch is the same. 

• The phenomenon leading to the transformation from anti-aligned to aligned periapses is 
clearly seen in the right-hand side panel, where the trajectory is seen crossing the origin of 
the plot. 

It is interesting to note that this picture has a counterpart in the study of periodic orbits 

of the 3-body problem. The study of symmetric periodic solutions shows the existence of two 
separated stable branches with aligned periapses; these two branches are tied with continuity by 
a branch of unstable periodic orbits (Hadjidemetriou, 2002). The unstable branch corresponds to 
the saddles shown in the fourth panel of fig. 6. In fig. 11 it would appear as a right shortcut on 
the horizontal axis tying the initial and final segments of stable solutions with aligned periapses. 

7 Conclusion 

The contents of this paper include with variable emphasis, the topics of a series of lectures whose 
main title was "Routes to Order: Capture into Resonance" . This was indeed the subject of the last 
section above. The study of this subject has however shown that differently from the restricted 
three-body problem, the capture into resonance drives the system immediately to stationary 
solutions known as "apsidal corotations" . The whole theory of these solutions was also included 
in the paper from the beginning, that is, from the formulation of the Hamiltonian equations of the 
planetary motions and the expansion of the disturbing function in the high-eccentricity planetary 
three-body problem. The secular theory of non-resonant systems was also given. Motions with 
aligned or anti-aligned periapses, resonant or not, resulting from non-conservative processes (tidal 
interactions with the disc) in the early phases of the system life, seem to be frequent in extra-solar 
planetary systems. 
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